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GLOBAL MAPPING OF EARTH-LIKE EXOPLANETS FROM SCATTERED LIGHT CURVES 
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ABSTRACT 

Scattered lights from terrestrial exoplanets provide val uable information about the planetary surface. 
Applying the surface reconstruction method proposed bv lFujii et alJ (|2010( ) to both diurnal and annual 
variations of the scattered light, we develop a reconstruction method of land distribution with both 
longitudinal and latitudinal resolutions. We find that one can recover a global map of an idealized 
Earth-like planet on the following assumptions: 1) cloudless, 2) a face-on circular orbit, 3) known 
surface types and their reflectance spectra 4) no atmospheric absorption, 5) known rotation rate 6) 
static map, and 7) no moon. Using the dependence of light curves on the planetary obliquity, we also 
show that the obliquity can be measured by adopting the minimization or the extended information 
criterion. We demonstrate a feasibility of our methodology by applying it to a multi-band photometry 
of a cloudless model Earth with future space missions such as the occulting ozone observatory (03). 
We conclude that future space missions can estimate both the surface distribution and the obliquity 
at least for cloudless Earth-like planets within 5 pc. 

Subject headings: astrobiology - Earth - scattering - techniques; photometric 



1. INTRODUCTION 

Recent progress in observational techniques has re- 
vealed various physical properties of exoplanets be- 
yond orbital parameters and planetary mass. Detec- 
tions of the atmospheric components have been reported 
for several systems using spectrosco py at the planetary 
transit and secondary eclipse (e.g. Charbonncau ct alj 
2002: Vidal-Madiar ct al.. 2003. 2004; Tinctti ct al. 2005 
Swain et al l200ail2009i) .' Interior compositio ns can be in- 



ferred from planetary mass and radius (e.g. iLeger et al.l 
120091 : ICharbonneau et al.ll2009[ ). Constructions of ther- 
mal maps of the planet ary atmosphere have been pro- 
posed by (IWilliams et al . 2006: Cowan & Agol 200$. A 
longitudinal thermal map o f HD 189733b has been con- 
structed by iKnutson et al.l (12007') based on the method 
proposed bv lCowan fc Agoll ()2008h . 

Nevertheless, an identification of planetary surface 
components still remains an ambitious challenge. One 
of the promising approaches is to use the scattered light 
of e xoplanets through the direct imaging observation 
(e.g.jSeager et al.ll2000l: iFord et al.ll200lHSudarskv et al.l 



120051 ). IFord et al.l (|2001| ) focused on the inhomogene- 
ity of the Earth surface which causes diurnal varia- 
tion of the scattered light. They computed the scat- 
tered light from a model Earth observed at a dis- 
tance of 10 pc and showed that time variations of the 
scattered light in different photometric bands highly 
depend on the geological and biological features on 
the planetary surface such as ocean, land, and even 
vegetation. More detailed characterizations (includ- 
ing spectroscopy) of the scattered light of the Earth 
and its time variations are discussed bo t h via Earth- 
shine observations (e.g. IWoolf et al.l 120021 : lArnold et al.l 
120021: iMontafies-Rodriguez et al.l 12006^ and simulations 



I^UU^ IMontaiies-Kodriguez et al.l izUUbl) and simul ations 
(|Tinetti et al.ll2006al |bl: lMontanes-Rodriguez et al.ll2006f ). 
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These studies have suggested a future possibility to in- 
vestigate the surface of Earth-like exoplanets by the scat- 
tered light curves. We note that such time variations of 
the scattered light are also appli cable to det e rmine the 
rotation period from as shown bv iPalle et al.l (|2008[ ). 

A variety of inversion techniques of the planetary 
surface from the scattered light curves have been pro- 
posed. Surprisingly, the first theoretical study of scat- 
tered light curves to make albedo maps has been per- 
formed at the beginning of the twentieth century al- 
though the author as su me asteroid and sate llites for the 
target (jRusselll [T906[) . ICowan et all (|2009[) performed 
principal component analysis (PGA) on multi-band pho- 
tometric data of the Earth observed by EPOXI (Ex- 
trasolar Planet Observation and Deep Impact Extended 
Investigation) mission, and extracted spectral features 
which roughly correspond to land and ocean. They 
also checked the time variation of these components and 
translated it to the longitudinal distri bution of these 
components based on the form ulation bv ICowan fc Agoll 
(|2008f) . lOaklev fc Cashl (|2009D paid attention to the gap 
of reflectivity between ocea n and land , and reproduced 
a longitudinal map of land. iFuiii et all (|2010f ) (hereafter 
FIO) have developed a methodology to estimate the areas 
of ocean, soil, vegetation, and snow from multi-band pho- 
tometry, and showed that the area of these components 
can be recovered from mock observations of a cloudless 
Earth. 

Since these authors focused on diurnal variations in 
mapping the surface, the resultant maps have only lon- 
gitudinal resolution with little information of the latitu- 
dinal distribution. One of the goals of the present paper 
is to develop a method to map inhomogcneous surfaces 
of exoplanets with both longitudinal and latitudinal res- 
olutions using both diurnal and annual variations of the 
scattered lights. 

We also consider the determination of the planetary 
obliquity from time variation of planetary light. The 
obliquity is an important property of Earth-like plan- 
ets with its strong implications for climate, habitabil- 
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ity ( e.g., iWilliams fc Kastind [l997t iWilliams fc PoUardl 
|2003[) and planetary formation. N-body simulations 
of the final stage of terrestrial planet formation indi- 
cated that th e distri b ution of the obliqu i ty C is isotropic 
(lAgiior et all [1999 !: Chambcrf^ [Ml IKokubo fc Idal 
120071 ) . iGaidos fc Wilhams (200l modeled the infrared 
light curves of exoplanets an d showed how the obli quity 
affects an annual variation. lOaklev fc CashI (|2009[ ) also 
pointed out a possibility to determine the planet's obliq- 
uity. In this paper, we demonstrate that the obliquity 
is estimated simultaneously with the global map of the 
planet by analyzing the scattered light curves over its 
orbital period. 

The rest of the paper is organized as follows. We first 
review the estimation method of the weighted area from 
multi-band photometry proposed by FIO, and describe 
our methodology to reconstruct the planetary surface 
and the obliquity measurement in §2. Assuming a fu- 
ture satellite mission for the direct imaging of Earth-like 
planets, we apply our methodology to mock observations 
based on real data of the scattering properties of the 
Earth in §3. Finally we summarize our results in §4. 

2. METHODS 

Let us briefly summarize the reconstruction method of 
the planetary surface from scattered light curves by FIO. 
Scattered light from the planetary surface is character- 
ized by the bidirectional reflectance distribution function 
(BRDF), /(z?o,'?i,'/';A)[str-i], where i?o, ^i, and ip are 
the incident zenith angle, the scattering zenith angle, and 
the angle between the incident and scattered light, re- 
spectively (Fig. [1] a). The BRDF is also a function of 
the wavelength A. Using the BRDF at (0, 6) fixed on 
the sphere ( Fig. [Tjb), the total scattered intensity at a 
given phase is provided by 

/(A) 



Under the above assumption, equation ([1} reduces to 



F*WRp / /(<^,0;A)cos^9o(0,^)cos^9l(0,^^)c^^], 

(1) 

where -F*(A) is the incident flux at wavelength A, Rp is a 
planetary radius, Syi is the surface area facing on the ob- 
server and illuminated by the host star (we call it the visi- 
ble and illuminated area), and dfi = sin 9d9d(j). The posi- 
tion on the surface (0, 0) fully specifies i?0: "(^i) and ip, and 
thus we write /(0, 6; A) = f{M<P, M^, ¥'(0,.^); A). 

In order to make a fitting model for the multi-band 
photometric data, FIG classified the planetary surface 
into four types: ocean, snow, soil and vegetation. Denot- 
ing the specific BRDF of each type k by fk{^o, '^i, 'p), 
one can write the local BRDF at position (0, 9) as a sum- 
mation of the specific BRDFs of the fc-th type surface 
weighted by local cover fractions m'-'^ -'((/), 6*): 



Nc 



/(0,0;A)==^A(,/),0;A)mW(0,0), 



(2) 



fc=i 



where Nc is the number of surface types and Nc = 4 in 
this case. Assuming that the scattering is isotropic (Lam- 
bertian), they approximate the BRDF by wavelength- 
dependence of albedo afc(A), 



afc(A) 



Nc 



I{X) = F4X)RlY, 



afe(A) 



k=l 

X f m^''\(l),e)cosM'P:&)cos-di{(j),0)dn. (4) 

5vi 

Assuming that they know spectra for each components, 
FIO developed a reconstruction method for the weighted 
area of the fc-th component. 



Ak = Ri 



w{4>, 0)m'^''\(j3, e)dvt 



S'vi 



_ cos'do{(f>,6) cos'di{(l),9) 



(5) 



where W{(j>,d) is the weight function. The above def- 
inition reduces equation ([J]) to a set of linear discrete 
equation: 



C 



Nc 



=:^afc(Af,)Afe, 



(3) 



Ceett / cosi?o('^,£')cosi?i((/),6l)drj, (6) 

where C is a factor that depends on the phase angle, 
i.e., the planctocentric angle between the host star and 
the observer, and Xb is the center wavelength of the 6-th 
band. Using a cloudless model Earth, they showed that 
the weighted area can be approximately recovered by 
solving equation ^ with an inequality condition Ak > 0, 
and the summation of Ak(t) is approximately constant 
over time, Ak(t) « const. « R^. It means that one 
can roughly estimate Rp as well by their method. There- 
fore, we use the fractional area Ak/Rp normalized by R^ 
in what follows. FIO have considered the spin rotation of 
a planet only and recovered the weighted area during one 
day of the planet and the weighted area can be converted 
to longitudinal maps by the i nversion method pr oposed 
bv lCowan fc"Agoil pOM ) and lCowan eTall (|2009[ ). 



2.1. An Inverse Problem of a Two-dimensional Map 

In the present paper, we consider the orbital motion 
to recover a two-dimensional world map of the planetary 
surface. We assume a model planet on a face-on circular 
orbit, which is the most promising case for the planet 
mapping. We also assume that orbital period is Port = 
365 [day] and spin rotational period is Pspin = 24 [hrs]. 
We denote the phases of the orbital motion and spin 
rotation by O and $ (Fig. [l] c), respectively. Our mock 
observation is assumed to be performed over one year 
from Q = Qs to 8 = Qs + 27r. There are two unknown 
parameters in this geometry: the obliquity and the 
initial orbital longitude 8s- We discuss the measurement 
of ( and 85 in section [2731 

For our conventional geometry, the weight function is 
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Figure 1. Schematic configurations of tiic planetary surface and system. Panel a) displays the definitions of the arguments of the BRDF 
f{'do,'dl,ip;X): the incident zenith angle, i9oi the zenith angle of scattering, and the relative azimuthal angle between the incident 
and scattered light, ip. Panel b) shows the spherical coordinate fixed on the planetary surface: the complementary latitude, 6, and the 
longitude, ip- The spin rotation axis is indicated by 6 = 0. The spherical coordinate of the vector from the planetary center to observer 
is denoted by ((/lobs i ^obs ) ■ Then, the spin rotation is described by ■!> = 27r — i^obs- Panel c) illustrates the coordinate system to specify 
the phase of the planet. The denotes the orbital longitude measured from the summer solstice (i.e. the angle between the incident ray 
and the projected vector of the spin rotation axis to the orbital plane). The obliquity f is the angle between the spin rotation axis and a 
normal vector of the orbital plane. In this paper, we assume that the line of sight is perpendicular to the orbital plane. 



given by 

W{(l), 61; 9; $; C) = ^ cos ?9o cos 
3 

= - [ cos 6 cos 9 sin C — sin 6 sin (0 + <I>) sin 9 

— sm 9 cos {(p + $) cos 9 cos C ] 

X [sin 6 cos {(j> + $) sin ( + cos 6 cos C] ■ (7) 

We provide the explicit form of the visible and illumi- 
nated area, Syi in Appendix A. Note that the weight 
function is normalized to unity by integrating over the 
sphere: 



J W^(0,0;9;$;C)dl^ = 1. 



(8) 

Figure [2] shows the behavior of the weight function 
W^(0, 61; 9; $; C) for cases of C = 0°, 45° and 90°. The 
weight function covers the region of 0° < < 90° -I- ( 
and 0° < (p < 360° as the planet rotates around the host 
star and its spin axis. As a result, the scattered light 
curves contain information about the land distribution 
up to 50 (1 + sinC)% of the total surface area in princi- 
ple. 

We consider how to deduce the local cover fraction 
m'*'^ (0, 9) if a set of the weighted area throughout the 
orbital and spin rotational periods is given. In order 
to solve this problem, we apply one of the techniques 
of tomography, the linear inverse problem, to the data. 
We assume that the planetary surface consists entirely of 
land and ocean and do not consider the other components 
nor further decompose them into clouds, vegetation or 
soil in what follows. Therefore we use the land cover 
fraction m(^, 9) for the surface classification {k = land). 
The surface at (0, 9) covered by land only (ocean only) is 
expressed by m{4),9) = 1 {raQ),9) ~ 0). We also denote 
the weighted area of land by A($(ii), 9(ii); C)- In §3, we 
also consider the components of vegetation and soil using 
mock photometric data. 

To begin with, let us discretize the observing time in 
N epochs and pixelize the planetary surface in M pix- 
els. The weighted area recovered at the i-th epoch ti is 



written as 

A(<i>(tO,9(tO;q 

= J2 w^(0,0;9(t,);$(iO;C)m(0,0)di7 

= E {m)^, f W{cP,9;e{U);<t>{U);C)dn,{9) 

nsvi)#0 

where Sj indicates the j-th pixel on the surface. The 
weighted cover fraction (m),;j is defined as 

_ i^Wi(P,9;e{t,);^uy,C)mi(f>,9)dn 
' J^Wi^,9;Q{U);^U);C)dn ' 

Here, let us define the pixel-averaged cover fraction 



TO, 



/^^ to(0, 9)dn 

J S i 



(11) 



Under the assumption that the weighted cover fraction 
approximately equals to the pixel-averaged cover frac- 
tion: 



m W TO 



3 I 



(12) 



we obtain 

A(<f>(t,),9(t,);C)=E^^ / w{cp,9;e{uy,^uyx)dn. 

(13) 

We introduce the data vector of weighted area a, the 
design matrix G and the model vector of cover fraction 
m as 

a = {a, = AifiU), eiUyO/cT^ U = 1, 2, N}, 
G = {G,,= / W{cl),9;e{t,y,^uyX)di}/a, 



I i = l,2,...,iV,j = l,2,...,Af}, 
m = {mj \j = 1,2,...,M}, 



(14) 
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Figure 2. Geometry dependence of the weight function VF{0, 6; 0; <3?; f ) in the visible and illuminated area 5vi- Annual variations 
are indicated by different thin lines. Solid, large-dashed, middle-dashed and small-dashed lines indicate = 0°, 90°, 180°, and 270°, 
respectively. Because the weight function solely depends on 0-1-$, in other words, longitude minus the phase of the spin rotation, we adopt 
(/) -|- to x-axes instead of showing diurnal variations explicitly. 

Thick lines indicate the boundary of the visible area, 9 = dv(<t>', Q in equation l)A3|l . Each panel shows a different obliquity, = 0° 
(panel a), 45° (panel b), and 90° (panel c). Inner and outer contours indicate VK(0, 6; 0; f ) = 0.6 and 0.3, respectively. 



where (Ti is the error of the weighted area. Then we can 
rewrite equation ()13p as 



a = G m. 



(15) 



The above equation can be solved for m by minimizing 

= Ifldata - Gm|^, (16) 

under inequahty constraints: 

0<mj<l, (17) 

where fldata is the observed weighted area with noise. We 
denote m that minimizes equation (|16p by most- 



2.2. Solving the Inverse Problem 

We are now in a position to solve equation ([TB| with 
the condition (fT7|) using an idealized data set. Since one 
can choose an arbitrary division for the surface modeling, 
we use the following pixelization: 

(l + sinC)A 



s,={{<P,6)\l 



Mg 



< COS < 



y\ ^ 1 27r 27r 
^^^A)— r;— , Yr(j0 - 1) < < —j<t>}^ 

Ma M, 



for j0 = 1,2,,, and jg = 1, 2, 



, Mg. 



(18) 



The above pixelization has a constant solid angle for a 
given C and the number of pixels does not change for 
different C's. We adopt 9 and Mg = 23 (see the 

right panel of Figure [3] for the C, = 90° case) and assume 
the numbers of epochs for diurnal variations -/V$ = 23 



and for annual variations Nq = 23. In Appendix B, 
we discuss the dependence of model degeneracies on the 
pixelization. 

We create the synthetic data of the weighted area us- 
ing 1° X 1° fixed land/water masks of the ISLSCP II 
{input map] left panel of Figure [3]). Pixel averaged land 
cover fraction of this data {reference map) is shown in 
the right panel of Figure [H We fix the obliquity Q = 90° 
and 0(to) = ©s = 0° and compute the weighted area 
^input('i', 0; C = 90°) from the input map (not from the 
reference map) . We compute the data Adata by adding a 
Gaussian noise Ng to Ainput('&, ©; C = 90°): 

Adata($, e) = Ainput(«', 6; C = 90°) + Ng. (19) 

We denote the standard deviation of the Gaussian ran- 
dom variable A'^g by cr = (N^)^/^. 

The linear inverse problem (eq. [15]), or the equiva- 
lent least square problem (eq. [H]), under the inequal- 
ity constraints (eq. [17| ) can be solved with Bounded 
Variable Least Squares Solve r (BVLS) developed by 
ILawson fc Hamol] ([1971 [1991 Q The BVLS is a gen- 
eralization of_yie_jiom^iegatiye least squar e (NNLS) de- 
scribed in ILawson &: HansonI (11974 |1995D and uses the 
QR decomposition (see also iMenkeJ 119891 ) . The BVLS 
works for both overdetermined and underdetermined (ill- 
condition) problems. Even for ill-condition problem, 
the BVLS provides one of the (non-unique) solutions 
(jLawson fc HansonI [l995l ). 

By solving the least square problem (eq. [TB]) with 

1 The original c ode of th e BVLS is available through NETLIB 
(jhttp:/ /www. netlib.org/lawson-hanson/indcx. html ) 
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Figure 3. The input map (left) and the referenc e map_( riKht). The left panel shows the land distribution of the Earth based on the data 
set of the ISLSCP II fixed land/water masks (ht tp:/ /islscp2.sesda.com| . Right panel indicates the pixel-averaged land cover fraction (Wij) 
computed by averaging the left panel over pixels. 



the BVLS algorithm, wc obtain the recovered map most ■ 
Figure |4] displays the two-dimensional images of the 
weighted area (left), the recovered maps (middle), and 
the prediction errors, [^data($(ii), 6(ii)) - Gmcstl/o- 
(right) for Cinput = 90°. In this figure, we regard 
^ and Os as fixed parameters and use the input val- 
ues C = Cnput = 90° and Qs = 6s,input = 0° for 
equation ([T4|). We consider four cases of the noise, 
a = 0, 0.01(Ainput),0.1(Ainput), and 0.3(Ainput), where 
(^input) is the average of Ainput(0, d;( = 90°) and is ap- 
proximately equal to 0.3 for the case of the Earth. As 
discussed in Appendix C, we confirmed that the esti- 
mated map is the unique solution under the inequality 
constraint. Comparing the recovered maps with the ref- 
erence map (right panel in Fig. [S]), we conclude that one 
can approximately recover the land distribution from the 
two-dimensional map of the weighted area. 

We also solve the inverse problem in the case of the 
obhquity of the Earth, Cinput = 23.45° (Figure [5|). Even 
for a small obliquity like the Earth, one can recover the 
land distribution of 50 [1 + sin (C = 23.45°)] « 70 % of 
the planetary surface. We conclude that our methodol- 
ogy can recover main features of the continents on the 
planetary surface even if the planet has small obliquity 
as is the case of the Earth. 

2.3. Measurement of Obliquity 

So far, we have assumed that the obliquity C and the 
initial phase Qs in the orbit are known a priori. Next 
we try to estimate the obliquity C by our method itself. 
Figure [6] indicates the obliquity dependence of images of 
the weighted We determine the best-fit value of C 

and 05 by minimizing distance: 

G(c,es) = {G,,(c,e5)= / w{c^,e;e{t,y,^t,yx)dn/a, 

J Sj 

|i = l,2,...,iV,j = l,2,...,M}, 

e(to) = es, (20) 

where m^g^" ^ is determined by minimizing with given 
C and 05 using the BVLS fitting. Figure [7] displays 
x2(C,es) for Cinput = pl^nd^inEut = 180°. We 
use the amoeba routine (jPress et al.lll992D to search for 
C and 05, which minimizes X^(C:0s)- We estimate er- 
rors of C and 0s by the bootstrap resampling because 



errors are originated not only from additional Gaussian 
noise but also from pixelization, linearization, and other 
systematics. After 10'^ iterations of the bootstrap, we 
obtain Cost = 46.3° ± 0.6° and 05,cst = 181.8° ± 0.8° 
for 10 % additional noise and Ccst = 50.3° ± 2.3° and 
05_o,t = 183.0° ± 2.3° for 30 % additional noise. The 
estimated values of both obliquity and initial phase in 
orbit remarkably agree well with the input values. 

3. APPLYING TO MOCK OBSERVATION OF MULTI-BAND 
PHOTOMETRY 

Now, we apply our methodology to a mock multi-band 
observation of the cloudless Earth as a more realistic 
demonstration. We simulate scattered light curves of a 
cloudless Earth in the same manner as the simulation 
in FIO. We use the solar spectrum as an incident flux 
on the model Earth. We assume an observer seeing the 
planet on a face-on and circular orbit at a distance of 
5 pc (Case A) or 10 pc (Case B) from it. The scat- 
tered light from land is computed using the BRDF of 
actual land surface assigned by MODIS data-set ( "snow- 
free gap-filled MODIS BRDF Model Parameters"). We 
adopt the data of April and ignore the seasonal variation 
of the BRDF at each point on the planetary surface. The 
scattered light from ocean is calcul ated with the BRDF 
model for wavy ocean described in iNakaiima fc Tanakal 
([1981 . In our model, the snow cover regions around the 
poles are replaced by ocean. We also include the effect 
of Rayleigh scattering by atmosphere with single scat- 
tering approximation between the atmosphere and the 
underlying surface. Neither the effect of clouds nor the 
molecular absorption is, however, incorporated. 

Our assumed observation parameters are listed in Ta- 
ble [TJ As an observing system, we assume a satellite 
telescope with 1.1 m aperture. This assumption comes 
from the basic architect ure proposed for th e occulting 
ozone observatory (03) (jKasdin et al.l [2010f ). which is 
a satellite mission for direct imaging of exoplanets in 
UV/optical/near-IR bands shading the light from the 
host star by a 30 m external occulter. The scattered light 
is computed at the central wavelength of each band and 
multiplied by its band width. We consider four bands 
centered at the wavelengths listed in Table [51 which cor- 
respond to the bands of the MODIS land BRDF data. 
We assume 0.1 /xm as the bandwidth. 

We divide the orbital phase into 23 (0^ = 27ri/23, i — 
0, 1, 22), and the spin rotational phase into 23 ($; = 
27ri/23, i — 0, 1, 22). We compute the scattered light 
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Figure 4. Two dimensional images of the weighted area A{4>, 9) (left panels), recovered maps (middle panels) for f = 90°, and prediction 
errors, [A^ata('J'(*i)i ~ Gmcst]/o" for different noises, a = 0, O.Ol(Ainput), 0.1(j4input), and 0.3(Ai,iput) from the top to the bottom 

panels. 
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Figure 5. Two dimensional images of tiie weiglitcd area A{4>,8) (left panels) and recovered maps (right panels) for f = 23.45° 
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Figure 6. Obliquity dependence of two dimensional images of the weighted area. 
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Figure 7. X^(Cj0s) for Cinput = 45° and ©s^input = 180° (cross) with 10 % additional noise. The estimated obliquity and initial orbital 
longitude are Cost = 46.3° ± 0.6° and es,ost =' 181.8° ± 0.8°, respectively. 

curves in 4 bands with the exposure time ioxp = 24/23 
hours throughout one year. Then, we stack data for 14 
days centered at each orbital phase &i, and fold the light 
curves according to its spin rotation period to obtain 
diurnal curves at the orbital phase, assuming that the 
spin rotational p eriod is already known through pcri- 
odogram analysis (jPalle et al.ll2008l ). Thus, the resultant 
light curves have 23 x 23 data points for each band. 

As for the noise, we consider the photon shot noise, the 
exzodiacal light, dark noise, and read noise. We ignore 
the leakage of the light from the host star. Thus, the 



signal to noise ratio S/N is expressed as 



S/N^^= ^ (21) 



where S, N^,, N^, and are the scattered light from the 
planet, the contribution from zodiacal light, dark noise, 
and read noise, respectively. Specifically, they are given 



Global Mapping of Exoplanets 9 



Table 1 

Observation parameters 



Synibol 


Quantity 


Value 


Unit 


texp 


exposure time 


24/23x3600 


sec 


Tl 


folded days 


14 




£) 


diameter of telescope aperture 


1.1 






sharpness 


0.0433 










SI vr''Sf^c 1 T\'\ Yol 


h 


end-to-end efficiency 


0.5 




V 


dark rate 


0.001 


counts/sec 


K 


read noise 


2 


\/ counts / read 


QE 


quantum efficiency 


0.91 






zodiacal light in magnitude 


23 


mag/square arcsec 


Z 


zodiacal light 


1 




Fo 


zero flux 


1.4 X lO"! (band 1) 


cts/cm-^/nm/s 






9.6 X 103(band 2) 


cts/cm'^/nm/s 






7.2 X 103(band 3) 


cts/ cni'^ / nm / s 






4.1 X 103(band 4) 


cts/cm^ /nm/s 



Table 2 

Band definition. 



band 


wavelength 


band width 




Afc [lira] 


A A [^m] 


1 


0.469 


0.1 


2 


0.555 


0.1 


3 


0.645 


0.1 


4 


0.8585 


0.1 



as 



1 

'2.5 _ TT 



Vtr. 



7V, = . 



(22) 
(23) 



(24) 



The parameters and their va lues we adopt are listed in 
Tablc[T](see e.g.. lSavranskv et al. 201Q). The zero flux Fq 
for calibration is calculated by linearly inte rpolating the 
spectr um of Vega given by Tables 1 and 2 in lColina et al.l 
(|1996[ ). and shown in Table [3] in this paper. In order to 
quantitatively assess each noise level, we list the ampli- 
tudes of these noises per epoch in Table [3] together with 
the total average of the signal S in the case of C = 90° . 
The averaged signal-to-noise per epoch is also listed in 
Table [3l We note that the averaged signal per epoch is 
comparable with the noise level. 

The obtained mock data l{\h) at each epoch 
(8(ti), $(ti)) is decomposed into the weighted area of 
four surface types, ( ocean, vegetation, soil, and snow 
) by solving equation ([6]) following FIO. The weighted 
area of land is computed by summing up that of soil, 
vegetation, and snow, that is, A(9, $) = ylsoii(0,$) -I- 
-4vegetation(S, $) + Asnow(6, $)• While FIO used nonlin- 
ear fitting in order to constrain to be positive, we 
use the BVLS linear solver instead without upper con- 
straints. The resultant A(8, $) as a function of Q and 
$ is displayed in the left panels of Figure El 

We now apply the reconstruction method described in 
§2 (see eqs [fl] [H] [17] ) to the set of obtained ^(6,$). 
For the errors Gi in equation (|14p . we first compute the 
variance of the weighted area by 100 times resampling of 



photon counts according to Poisson statistics, and we call 
such variance as erf j.^^. The Ci^rcs, however, turned out 
not to relevant for Gi in equation (jl4p because it often 
results in zero or very small errors when Ai = due to the 
boundary condition {Ai > 0). Therefore, we substitute 
an average value of errors a = J2i=i <^i,res/N to equation 
instead. For the case that all estimated errors of Ak 
are positive, the results with cr^ = ai^^es do not change 
significantly. In the fitting, we assume that the local 
cover fraction of each surface type remains unchanged 
and thus ignore the seasonal variations of vegetation, soil, 
and snow. The center and right rows of Figure [8] display 
recovered maps, and prediction errors for Cinput = 90° of 
the case A and B. The major features of the surface can 
be seen in the resultant maps. The averaged prediction 
errors | Adata($(i»), e(t,)) - Gm^st\/N are 0.53 (Adata) 
(case A) and 0.85 (Adata) (case B) , respectively. 

The land and ocean distributions have been consid- 
ered so far. Next we create color composition maps of 
vegetation and soil on the land and ocean recovered map 
(Figure [9]). Adopting the same method of the land recov- 
ery to the weighted areas of soil {Asoi\{ti)) and vegetation 

(^vegetation 

(U)), we derive the soil and vegetation distri- 
butions and blend them by yellow (soil) and green (veg- 
etation) over Figure m in other words, the land (white) 
and ocean (blue) distribution. Even using the weighted 
areas of soil and vegetation separately for the recovery, 
the resulting vegetation and soil maps almost distribute 
over the land that recovered with the summation of the 
weighted areas (Aland (ii) = AsoniU) + Avegetation(ti) + 
AsnowiU))- In this figure, snow is ignored because the 
contribution of it is significantly small. The concentra- 
tions of vegetation at South America and soil at middle 
Africa seen in the case A indeed correspond to the Ama- 
zon forest and the Sahara desert, respectively. 

Now, we try to measure the obliquity for the mock ob- 
servation. Figures [TUl [TTl [TH and [T3] display the light 
curves of 4 bands with 6s input = 180°, corresponding 
to difl^erent obliquities, Cinput = 90°, 60°, 45°, and 30°, 
respectively. The noise assuming in these figures are 
computed in the "Case A" . Performing the fitting de- 
scribed in §2.31 we estimate the best-fit obliquities and 
the initial orbital longitude. The predicted curves with 
the best- fit values are shown by solid lines in Figures [TUl- 
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Table 3 

Assumed signal and noises in 4 bands. 



band b 



1 



averaged signal S (case A) [10^ cts/epoch] 
averaged signal S (case B) [10^ cts/epoch] 
exozodi Nz [10^ cts/epoch] 
dark noise [10^ cts/epoch] 
read noise [10^ cts/epoch] 
averaged S/N (case A) [1/epoch] * 
averaged S/N (case B) [1/epoch] * 



40.9 
10.2 
12.4 



3.3 
0.83 



40.3 36.9 
10.1 9.2 
9.2 7.0 

1.2 

1.3 
3.6 3.7 
0.93 0.94 



40.9 
10.2 
4.0 



4.9 
1.3 



* The averaged S/N are defined by S/y'S + N^ + Na + Nr 
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Prediction Error in units of a 





^^^^^^^^^+— 









0.2 0.4 0.6 0.8 


H 1 

1 






0.2 



0.4 



0.6 



0.8 





0.2 0.4 0.6 0.8 
(* + &)/ 2 n 

Weighted Area 



Case A (5 pc), <-=90.00°, Band 1-4 
Land cover fraction 



0.2 0.4 0.6 0.8 1 

(* + e) / 2 77 

Prediction Error in units of u 




—I — I — I — I I — I ^ i i — 



0.2 0.4 0.6 0.8 

(* + 0)/ 2 TT 



Case B (10 pc), <-=90.00°, Band 1-4 



0.2 0.4 0.6 O.f 
(* + 8) / 2 77 



Figure 8. Simulated weighted area from mock photometric observation (left), its recovered map (middle), and prediction error (right). 
Top and bottom panels indicates the case A and B, respectively. We note that the snow component such as the South Pole in the input 
data is replaced to the ocean. 



[T51 The input and the best-fit vahics and the reduced 
are fisted in Table ID We compute the for the land dis- 
tribution (Eq. [20] )■ Therefore, the degree of freedom is 
23 X 23-23 X 9 = 322. It is likely that the inequality con- 
straints and the other systematics such as the Lambert 
assumption induce the relatively high reduced x^ values 
(~ 2). One can see both diurnal and annual variations 
of light curves in these figures. The larger variations are 
seen in band 4 ( 0.8585 yU m) due to large albedos of 
the land component (soil and vegetation; see Figure 7 in 
FIO). There is "pinching" of the light curve at 85 = 0° 
(i ~ 264) for Cinput = 7r/4 because the maximum of the 
weight function is located at the planet's pole. In this 
period, the physical position of the maximum stays at 
the pole for the daily motion. As a result, there are little 
variations in the period. However the "pinching" seen at 
65 - 180° (i - 0) for Cinput = 90° has a different ori- 
gin. In this period, the weight function moves around 



the Southern Hemisphere, that is, ocean mainly. 

The maps for the four input obliquities are dis- 
played in Figure [TJ] and the best-fit obliquities and ini- 
tial longitudes are listed in Table |4l For the case A, 
the fitting provides the estimated values close to the 
input ones within 10°. Figure fTSl shows the difference 
of the predicted light curves with the best-fit obliquity 
Cost = 44.4 and that shifted ±3.3° (~ la level estimated 
by the bootstrap resampling), C = 47.7° and 41.1° . 

The results for the case B are also listed in Table U) 
Figure [16] displays the map for the case B. The es- 
timated values for Cinput = 30° and 60° tend to deviate 
from the input one. As shown in the bottom right panel, 
the result for Cinput = 30° are significantly biased. It is 
likely that the systematics affect the obliquity estimates 
for this noise level. We also perform the other method 
called the extended information criterion in Appendix D. 
Although the tendency that the estimated value is biased 
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Case A (5 pc), <-=90.00°, Band 1-4 Case B (10 pc), ^=90.00°, Band 1-4 




CLASSIFICATION BASED ON IGBP 

Figure 9. Recovered color maps of soil and vegetation from the mock observations. Soil and vegetation are blended by yellow and green 
on the land recovered maps showrn in the middle panels of Figure [8] Top left and right panels indicate the case A and case B, respectively. 
Snow is replaced to ocean in this simulation. Bottom panel is a classification map of ocean (blue), snow (pink), vegetation (green), and 
soil (yellow) based on the IGBP classification. Our classification of soil and vegetation is listed in Table 2 of FIO. 




Figure 10. The mock light curves in units of reflectivity for Cinput = 90°. Each panel indicates a different band: top left, top right, bottom 
left, and bottom right panels corresponds to bands 1 (0.469/^ m), 2 (0.555/.t m), 3 (0.645^ m) and 4 (0.8585^ m). Dashed vertical lines 
divide each epoch ■!>, stack in 14 days observations. There are 23 data points in one epoch and 23 epochs in whole data set. Therefore, 
there are 23 = 529 data points in whole data set (i = 1, 2, , , 529). Wider versions of two shaded regions labeled in a) and b) are inserted 
in the mini panels. The predicted curves with the best-fit obliquities (^est = 89.8° and initial longitude ©Scst = 183.8 are drawn by solid 
lines. The observational noises are computed on the assumption of case A. 

to C = 90° vanishes, the best-fit value does not improve significantly. 
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Figure 11. Same as Figure [TO] for Cinput = 60°. 
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Figure 12. Same as Figure [To] for Cinput = 45°. 

In short summary, the planetary obliquity can be es- 
timated well for the case A, while, for the case B that 
has larger noises, it is marginal to derive reliable values 
of the estimated obliquity by our method. 

4. SUMMARY AND DISCUSSION 

We have developed the reconstruction method of the 
two-dimensional planetary surface via diurnal and annual 




100 200 300 400 500 

i 




100 200 300 400 500 




variation of the scattered light. Applying the method to 
the mock photometric data, we have demonstrated that 
our method works for the mock Earth model, while this 
model has a lot of simplifying assumptions as follows: 1) 
cloudless, 2) a face-on circular orbit, 3) known reflectance 
spectra 4) no atmospheric absorption, 5) known rotation 
rate 6) static map, and 7) no moon. We also found that 



Global Mapping of Exoplancts 



13 




Table 4 

Estimated obliquity and the initial orbital longitude by minimization and the EIC. 









the fit 


the fit 






EIC 


EIC 


Case 


Cinput 


©S.input 


Ccst 


©S.ost 




xVdof 


Ccst 


6s,est 


A 


90° 


180° 


89.8° ± 1.2° 


183.8 ± 2.2° 


650.3 


2.02 


89.0° 


184.5° 


A 


60° 


180° 


64.7° ± 4.8° 


175.3 ±3.5° 


661.2 


2.05 


65.3° 


176.5° 


A 


45° 


180° 


44.4° ± 3.3° 


179.8 ±4.1° 


731.1 


2.27 


44.5° 


177.5° 


A 


30° 


180° 


32.5° ± 3.0° 


188.6 ± 5.6° 


677.5 


2.10 


30.7° 


190.5° 


B 


90° 


180° 


89.8° ±3.5° 


190.8° ± 17.5° 


596.1 


1.85 


88.1° 


187.5° 


B 


60° 


180° 


79.7° ± 8.1° 


171.8° ± 24.8° 


568.5 


1.77 


48.5° 


155.6° 


B 


45° 


180° 


38.2° ± 21.6° 


231.8° ± 43.6° 


609.5 


1.89 


26.7° 


132.6° 


B 


30° 


180° 


89.2° ± 17.1° 


210.3° ± 32.7° 


714.4 


2.23 


7.9° 


138.6° 



the planetary obliquity can be estimated by this method. 
With our method, future sat ellite missions such as the 
occulting ozone observatory (iKasdin et al.ll2010f ) might 
provide "a global map" of Earth-like exoplanets. While 
only the terrestrial planets have been considered in this 
paper, our method might be applicable to any planets 
with an inhomogeneous surface, including Jupiter-like 
exoplanets. 

In this paper, we ignored the effect of clouds and ex- 
pected that clouds affect the estimation as like a statis- 
tical noise because of relatively short time va r iation of 
clouds. The PGA performed by iGowan et all (|2009[ ) is 
one of promising approaches because they could separate 
the land and ocean compositions even though they used 
the EPOXI data that contains the cloud effect. The ef- 
fect of clouds is discussed elsewhere (Fujii et al. in prepa- 
ration) .We also assumed a face-on circular orbit in this 



paper. This assumption might be too severe for practical 
applications. Wc will generalize our method in the next 
paper. 
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ture. Sports, Science and Technology (Nos. 20-10466 and 
22-5467), and by the JSPS Gorc-to-Corc Program "Inter- 
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Figure 14. The Q and @s dependence of X^(Ci0s) of the mock observational data set (case A) for ©s, input = 180° and the different 
input obUquities, Cinput = 90° (top left) , Cinput = 60° (top right) , Cinput = 45° (bottom left) , Cinput = 30° (bottom right). The input 
values are indicated by cross. The estimated obliquity and initial orbital longitude arc listed in Table |4] 



APPENDIX 
A. VISIBLE AND ILLUMINATED AREA 
The visible and illuminated area is expressed as 

r {(0, e)\Q < e < min[0v{<P; c), 0iA<l>: ©; C)]} 
„ _ I for < e < 7r/2 or 3n/2 < 9 < tt , 

~ ] {(0, e)|07.-(<^; e; $; C) < e < ^^(0; C)} 
[ for 7r/2 < e < 37r/2. 

(Al) 

The boundary lines of visible area 6'v/(0; $; C,) and illuminated area 9i^+{(j); Q; $; (), and 9i-{(j); 9; $; C,) are obtained 
by solving the following equation. 

W^(0,0;e;$;C) 
= (sin cos (0 + $) sin C, + cos 6 cos 

X (cos 9 cos sin C, — sin 6* sin (0 + $) sin Q — sin 6* cos (</> + $) cos cos Q 

= 0. (A2) 
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Figure 15. The difTercnco between the light curve with the best-fit values (Ccst = 44.4°, ©s est = 179.8°), with a 3.3° shifted obliquity 
((J = 47.7°, Bs = 179.8°; solid curves) and with a -3.3° shifted obhquity (C = 41.1°, 65 = 179.8°; dashed curves) 



The explicit equations of 0^(0; $; C), ()i,+ {4>^ C)i and 9i^^{(p; 9; $; C) are expressed as 

cos {(f) + $) sinC 



6ly((/i;$;C)==cos" 



cos2 ( -f cos2 ((^ -f (()) sin^ C, 



9/,+ (0; 8; $; C) - cos-i (7K<^ + $; 9; 0) 
9,,_ (0; 9; $; C) = cos-^ (-77(0 + $; 9; C)), 



(A3) 

(A4) 
(A5) 



where 



,7(^;9;C)^ 



cos C COS 9 cos (j) + sin 9 sin c 



\J cos^ C cos^ 9 cos^ + 2 cos C, sin 9 cos 9 sin cos + sin^ C cos^ 9 + sin^ 9 sin^ 



B. PIXELIZATION OF PLANETARY SURFACE 

We note here that the design matrix alone does not fully specify mcst without any constraints. We perform singular 
value decomposition (SVD) of the design matrix: 

G = UkV^, (Bl) 
where U and V are N x N and M x M unitary matrices, and A is a N x AI matrix: 



A = 



A^ 



A = diag{li, ..Jm), 



(B2) 
(B3) 



with Ij being the singular value of G in descending order. Wc denote the j-th row of V by the vector Vj. 

In general, there are A'nuii zero components in the singular values, (Im-n^^u+i = Im-n„^u+'^ = ^ 0)- Corre- 
sponding V vectors are called null vectors. Linear combination of the null vectors do not affect the data vector, 



G J2 /3"""'''v"""''=0, 



null.s — 



(B4) 

where is an arbitrary coefficient. Equation (|B4p indicates that the predicted model has a freedom to add any 

linear combination of the null vectors. 
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Solid curve in the left panel of Figure [17] displays the singular values in descending order (this plot is referred to 
as spectrum of data kernel) in the case of = = Nq = 23 and Ng ~ 9. Due to the truncation errors, the 
singular values obtained by numerical algorithms of the SVD do not vanish exactly. However, only three of them 
{j — M ~ 2, M — 1, and M) have significantly smaller values than other diagonal elements. We regard the elements 
smaller than 10"'^ as zero and the corresponding row vectors of V as null vectors. 

Although the presence of null vectors does not directly lead to the indcfiniteness if one adopts the inequality 
constraints fea.|17p. smaller number of null vectors makes easier to interpret an uncertainty of the recovered map. 
Therefore it is useful to know how the number of null vectors depends on pixelization. Then one can decrease the 
number of the null vectors by an appropriate choice of pixelization. For simplicity, we fix = iV$ = A'e and estimate 
the number of the null vectors by changing Mg and M^. The left panel in Figure [T71 shows the spectra of the singular 
value in descending order for three cases; = 22 (dotted), = 23 (solid) and AI^ = 24 (dashed) with Mg = 9. 
The case of = 23 has a steeper spectrum than the others. It means that the choice of = 23 makes easier to 
identify null vectors. The right panel of Figure [17] indicates the number of null vectors for different sets of and 
Mg, where we define vectors with singular value below 10~^ as null vectors. There is a general tendency that the odd 
number bin of AI^ has lesser iVnuii- Although the selection of smaller values of Alg or generally leads to smaller 
number of the null vectors, it degrades the resolution of the map. As a compromise, we decided to adopt = 23 
and Mg = 9 for our fiducial values of the model in this paper. Figure [T5I indicates the (unit) null vectors for = 23 
and Mg = 9, u™!''! = vnj^2, t?"""'-^ = vm-i, and — where Vj is the j-th row of V. These three null vectors 

for Mff, = 23 and Mg = 9 affect the large-scale structure of the map but do not change the small-scale distribution like 
the continents of the present Earth. 
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Figure 17. Singular values of the design matrix G in descending order (left panel). Dotted, solid, and dashed curves correspond to 
AI^ = 22, AI^ = 23 and = 24, respectively. All curves assume Mg = 9. Right panel shows the number of vectors with singular values 
below 10^* as a function of Mg and AI^ = 7V$ = A^e- 






Unit null vector 1 



Unit null vector 2 



Unit null vector 3 



Figure 18. Three unit null vectors of our fiducial model A/^ = = N@ = 23 and AJg = 9. These maps are drawn using Hammer- AitofT 
projection. 



C. UNIQUENESS OF THE RECOVERED MAP 

In this appendix, we examine whether there is the freedom for adding any components of null vectors to the recovered 
map derived by the BVLS. Because of the boundary condition for mcst (eq. |17p. the estimated values of m^stj are 
classified into 3 types — (a) rricstj ~ 1 (on the upper boundary), (b) rricstj = (on the lower boundary), and (c) 
0< mestj < 1 (in between). We focus on types (a) and (b) to consider the uniqueness of our recovered map. We 
denote the index j of type (a) by j— and that of type (b) by j+. In order not to violate the boundary condition, the 
linear combination of null vectors should satisfy the following constraints: 

^ ^„un,fc^nun,fe<o for any (CI) 
fe=i 

3 

^null,fe^nun,fe > Q ^^^y .^^ ^^2) 

fe=l 

These inequalities can be solved in terms of using the fact that Vj is positive (see Fig. [18]): 

null, 2 null, 3 

^nulU < _!j_^null,2 _ !ij_^null,3 (^3) 
— null,l^ null.l^ ' ^ ^ 

^j- 

null, 2 null, 3 

^null,l > ^nun,2 _ !l+ ^nulL3^ (^4) 

" — nulLl ^ null.l ^ 

Then, the following inequality is required so that /3"""'i exists, 

null, 2 null, 3 null, 2 null, 3 

^i- onull,2 I ''^3- onull,3 ^ onull,2 , + 



/?-11^2 ^Z^^null,3 < ^+ ^uun,2 ^+ ^uun,3_ (^5) 
11,1'^ null.l — null,!'^ null,!'^ ^ 



nuli,_ 



This equation can be reduced to the following expression: 

P(,+,,_) • b > (C6) 
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where p(j_). j_) and b are 2-diniensional vectors: 



(null, 2 null, 2 null, 3 null, 3 \ 

null.l null,!' null,l null,l I ' ' 

^^^^null,2^^null,3^ (C8) 

Now let us denote the arguments of p and b by 9a{i+,j—) and Op, respectively, i.e., 

PO+J-) = |PO+J-)l (cos6'„(j+,j-),sin6'„(j+, j-)), (C9) 
b=|b| (cos6'^,sin6l/3). (CIO) 
If |b| 7^ 0, the solution of equation (jC6[) should satisfy 

Oo.{]+,:i-)<ep<e^{j+,j-). (cii) 

We calcul ate t he above constraints for any combination of ), and found that there are no On that satisfies 

equation (jCll() for all c ombinations of )■ Thus, the only solution allowed is _ ^nuil,3 _ From 

equations jCH and JCD), it follows that =0. As a resuh, we confirmed that there is no room for adding extra 

components of null vectors to the recovered map. 

D. ANALYSIS WITH THE EXTENDED INFORMATION CRITERION 

In this appendix, w e try to perform the model selection of Cost and 0s,cst using the extended information criterion 
(EIC; e.g. lEfronlHOSl iKonishi fc Kitagawalll996l: UshTguro et al.lll997l) based on the KuUback-Leibler divergence. Al- 
though the model selection for Q and 0s based on the minimization is very convenient, the result might be biased 
for large noises compared with signal. The EIC is a bootstrap-based extension of Akaike Information Criterion defined 

as 

AT 

S/C= -2^1og^(a|mest,Ccst,es,ost) + 26B, (Dl) 

1=1 

where J^(a|most, Cost, 0s,cst) indicates the likelihood function . The best m odel or parameter set are given by minimizing 
the EIC. The bootstrap estimate of the bias bs is given by (llshiguro et al.„1997i) : 

^ JVeic 

bs = ^r^ V [logJP(a*''''|m*;^ Cost'', ©sL) -log-^(a*'''l"^cst, Cost, es,cst) 

JVEIC 



fc=l 

+ log J'(a|mcst, Cost, 05,ost) - log J"(a|m*,t , Cst , ©s.ost)]' (D2) 

"ost , Cost , and 65 indicate the fc-th bootstrap rcsample, most. Cost, and Os^est, and iVsic is the number 
of trials of the bootstrap. Assuming a Gaussian distribution for a with (7 = 1, we derive 

EIC^\a- G(C, es)mo,t I' + N log 2^ + bg 



Weic 

2iVEic' 



^'B^-^n^ E[|"*''-G(C*'''©5V:;t -|a*'^-G(C,es)most| 



fe=i 

+ \a- G(C, e5)most|' - |a - G(C*^^ )m:;^ |']. (D3) 

We adopt A^eic = 500 for bootstrap resampling. Figure [19] display the dependence of the EIC on Cost and 6s,cst- 
The minimum EIC estimate listed in Table |4] is derived by a stepwise search at one degree interval. 

The minimum EIC estimate listed in Table |3] is derived by a stepwise search at one degree interval. For the case A, 
the EIC provides the almost same values as listed in Table H) Figure [19] displays the EIC map for the case B. Although 
the estimated obliquities for Cinput = 30°, 45° and 60° by the EIC do not agree well with the input value (Table 4), 
the bias of Cinput = 60° and 30° seems to be improved compared to that from fitting. 



Global Mapping of Exoplancts 
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Case B EIC Case B EIC 




20 40 60 80 20 40 60 

<■ [deg] <- [deg] 



Figure 19. The C, and 0^ dependence of the EIC of the mock observational data set (case B) for 05, input = 180° and the different input 
obhquitics, ^input — 90° (top left) , C input — 60° (top right) , ^input — 45° (bottom left) , ^input — 30° (bottom right). The input values 
arc indicated by cross. 
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